home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Diamond Collection / The Diamond Collection (Software Vault)(Digital Impact).ISO / cdr44 / newmat08.zip / NEWMAT3.CPP < prev    next >
C/C++ Source or Header  |  1995-01-16  |  16KB  |  569 lines
  1. //$$ newmat3.cpp        Matrix get and restore rows and columns
  2.  
  3. // Copyright (C) 1991,2,3,4: R B Davies
  4.  
  5. #include "include.h"
  6.  
  7. #include "newmat.h"
  8. #include "newmatrc.h"
  9.  
  10. //#define REPORT { static ExeCounter ExeCount(__LINE__,3); ++ExeCount; }
  11.  
  12. #define REPORT {}
  13.  
  14. //#define MONITOR(what,storage,store) \
  15. //   { cout << what << " " << storage << " at " << (long)store << "\n"; }
  16.  
  17. #define MONITOR(what,store,storage) {}
  18.  
  19.  
  20. // Control bits codes for GetRow, GetCol, RestoreRow, RestoreCol
  21. //
  22. // LoadOnEntry:
  23. //    Load data into MatrixRow or Col dummy array under GetRow or GetCol
  24. // StoreOnExit:
  25. //    Restore data to original matrix under RestoreRow or RestoreCol
  26. // IsACopy:
  27. //    Set by GetRow/Col: MatrixRow or Col array is a copy
  28. // DirectPart:
  29. //    Load or restore only part directly stored; must be set with StoreOnExit
  30. //    Still have decide how to handle this with symmetric
  31. // StoreHere:
  32. //    used in columns only - store data at supplied storage address, adjusted
  33. //    for skip; used for GetCol, NextCol & RestoreCol. No need to fill out
  34. //    zeros.
  35.  
  36. // For symmetric matrices, treat columns as rows unless StoreHere is set;
  37. // then stick to columns as this will give better performance for doing
  38. // inverses
  39.  
  40.  
  41. // These will work as a default
  42. // but need to override NextRow for efficiency
  43.  
  44. // Assume pointer arithmetic works for pointers out of range - not strict C++.
  45.  
  46.  
  47. void GeneralMatrix::NextRow(MatrixRowCol& mrc)
  48. {
  49.     REPORT
  50.     if (+(mrc.cw*StoreOnExit)) { REPORT this->RestoreRow(mrc); }
  51.     if (+(mrc.cw*IsACopy))
  52.     {
  53.         REPORT
  54.         Real* s = mrc.store + mrc.skip;
  55.         MONITOR_REAL_DELETE("Free   (NextRow)",mrc.storage,s)
  56. #ifdef Version21
  57.         delete [] s;
  58. #else
  59.         delete [mrc.storage] s;
  60. #endif
  61.     }
  62.     mrc.rowcol++;
  63.     if (mrc.rowcol<nrows) { REPORT this->GetRow(mrc); }
  64.     else { REPORT mrc.cw -= (StoreOnExit+IsACopy); }
  65. }
  66.  
  67. void GeneralMatrix::NextCol(MatrixRowCol& mrc)
  68. {
  69.     REPORT                                                // 423
  70.     if (+(mrc.cw*StoreOnExit)) { REPORT this->RestoreCol(mrc); }
  71.     int t1 = +(mrc.cw*IsACopy); int t2 = !(mrc.cw*StoreHere);
  72.     if ( t1 && t2 )
  73.     {
  74.         REPORT                                             // not accessed
  75.         Real* s = mrc.store + mrc.skip;
  76.         MONITOR_REAL_DELETE("Free   (NextCol)",mrc.storage,s)
  77. #ifdef Version21
  78.         delete [] s;
  79. #else
  80.         delete [mrc.storage] s;
  81. #endif
  82.     }
  83.     mrc.rowcol++;
  84.     if (mrc.rowcol<ncols) { REPORT this->GetCol(mrc); }
  85.     else { REPORT mrc.cw -= (StoreOnExit+IsACopy); }
  86. }
  87.  
  88.  
  89. // routines for matrix
  90.  
  91. void Matrix::GetRow(MatrixRowCol& mrc)
  92. {
  93.     REPORT
  94.     mrc.skip=0; mrc.cw-=IsACopy; mrc.storage=ncols; mrc.length=ncols;
  95.     mrc.store=store+mrc.rowcol*ncols;
  96. }
  97.  
  98.  
  99. void Matrix::GetCol(MatrixRowCol& mrc)
  100. {
  101.     REPORT
  102.     mrc.skip=0; mrc.storage=nrows; mrc.length=nrows;
  103.         int t1 = !(mrc.cw*StoreHere);
  104.     if ( ncols==1 && t1 )
  105.         { REPORT mrc.cw-=IsACopy; mrc.store=store; }
  106.     else
  107.     {
  108.         mrc.cw+=IsACopy; Real* ColCopy;
  109.         if ( !(mrc.cw*StoreHere) )
  110.         {
  111.             REPORT
  112.             ColCopy = new Real [nrows]; MatrixErrorNoSpace(ColCopy);
  113.             MONITOR_REAL_NEW("Make (MatGetCol)",nrows,ColCopy)
  114.             mrc.store = ColCopy;
  115.         }
  116.         else { REPORT ColCopy = mrc.store; }
  117.         if (+(mrc.cw*LoadOnEntry))
  118.         {
  119.             REPORT
  120.             Real* Mstore = store+mrc.rowcol; int i=nrows;
  121.             while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols; }
  122.         }
  123.     }
  124. }
  125.  
  126. void Matrix::RestoreCol(MatrixRowCol& mrc)
  127. {
  128. //  if (mrc.cw*StoreOnExit)
  129.     REPORT                                   // 429
  130.     if (+(mrc.cw*IsACopy))
  131.     {
  132.         REPORT                                // 426
  133.         Real* Mstore = store+mrc.rowcol; int i=nrows; Real* Cstore = mrc.store;
  134.         while (i--) { *Mstore = *Cstore++; Mstore+=ncols; }
  135.     }
  136. }
  137.  
  138. void Matrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrMat(); }  // 1808
  139.  
  140. void Matrix::NextCol(MatrixRowCol& mrc)
  141. {
  142.    REPORT                                        // 632
  143.    if (+(mrc.cw*StoreOnExit)) { REPORT RestoreCol(mrc); }
  144.    mrc.rowcol++;
  145.    if (mrc.rowcol<ncols)
  146.    {
  147.       if (+(mrc.cw*LoadOnEntry))
  148.       {
  149.      REPORT
  150.          Real* ColCopy = mrc.store;
  151.          Real* Mstore = store+mrc.rowcol; int i=nrows;
  152.          while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols; }
  153.       }
  154.    }
  155.    else { REPORT mrc.cw -= StoreOnExit; }
  156. }
  157.  
  158. // routines for diagonal matrix
  159.  
  160. void DiagonalMatrix::GetRow(MatrixRowCol& mrc)
  161. {
  162.    REPORT
  163.    mrc.skip=mrc.rowcol; mrc.cw-=IsACopy; mrc.storage=1; mrc.store=store;
  164.    mrc.length=ncols;
  165. }
  166.  
  167. void DiagonalMatrix::GetCol(MatrixRowCol& mrc)
  168. {
  169.    REPORT 
  170.    mrc.skip=mrc.rowcol; mrc.storage=1; mrc.length=nrows;
  171.    if (+(mrc.cw*StoreHere))
  172.       { REPORT *(mrc.store+mrc.rowcol)=*(store+mrc.rowcol); mrc.cw+=IsACopy; }
  173.    else { REPORT mrc.store = store; mrc.cw-=IsACopy; }     // not accessed
  174. }
  175.  
  176. void DiagonalMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrDiag(); }
  177.                               // 800
  178. void DiagonalMatrix::NextCol(MatrixRowCol& mrc)
  179. {
  180.    REPORT
  181.    if (+(mrc.cw*StoreHere))
  182.    {
  183.       if (+(mrc.cw*StoreOnExit))
  184.          { REPORT *(store+mrc.rowcol)=*(mrc.store+mrc.rowcol); }
  185.       mrc.IncrDiag();
  186.       int t1 = +(mrc.cw*LoadOnEntry);
  187.       if (t1 && mrc.rowcol < ncols)
  188.          { REPORT *(mrc.store+mrc.rowcol)=*(store+mrc.rowcol); }
  189.    }
  190.    else { REPORT mrc.IncrDiag(); }                     // not accessed
  191. }
  192.  
  193. // routines for upper triangular matrix
  194.  
  195. void UpperTriangularMatrix::GetRow(MatrixRowCol& mrc)
  196. {
  197.    REPORT
  198.    int row = mrc.rowcol; mrc.skip=row; mrc.cw-=IsACopy; mrc.length=ncols;
  199.    mrc.storage=ncols-row; mrc.store=store+(row*(2*ncols-row-1))/2;
  200. }
  201.  
  202.  
  203. void UpperTriangularMatrix::GetCol(MatrixRowCol& mrc)
  204. {
  205.    REPORT
  206.    mrc.skip=0; mrc.cw+=IsACopy; int i=mrc.rowcol+1; mrc.storage=i;
  207.    mrc.length=nrows;
  208.    Real* ColCopy;
  209.    if ( !(mrc.cw*StoreHere) )
  210.    {
  211.       REPORT                                              // not accessed
  212.       ColCopy = new Real [i]; MatrixErrorNoSpace(ColCopy);
  213.       MONITOR_REAL_NEW("Make (UT GetCol)",i,ColCopy) 
  214.       mrc.store = ColCopy;
  215.    }
  216.    else { REPORT ColCopy = mrc.store; }
  217.    if (+(mrc.cw*LoadOnEntry))
  218.    {
  219.       REPORT
  220.       Real* Mstore = store+mrc.rowcol; int j = ncols;
  221.       while (i--) { *ColCopy++ = *Mstore; Mstore += --j; }
  222.    }
  223. }
  224.  
  225. void UpperTriangularMatrix::RestoreCol(MatrixRowCol& mrc)
  226. {
  227. //  if (mrc.cw*StoreOnExit)
  228.   {
  229.      REPORT
  230.      Real* Mstore = store+mrc.rowcol; int i=mrc.rowcol+1; int j = ncols;
  231.      Real* Cstore = mrc.store;
  232.      while (i--) { *Mstore = *Cstore++; Mstore += --j; }
  233.   }
  234. }
  235.  
  236. void UpperTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrUT(); }
  237.                               // 722
  238.  
  239. // routines for lower triangular matrix
  240.  
  241. void LowerTriangularMatrix::GetRow(MatrixRowCol& mrc)
  242. {
  243.    REPORT
  244.    int row=mrc.rowcol; mrc.skip=0; mrc.cw-=IsACopy; mrc.storage=row+1;
  245.    mrc.length=ncols; mrc.store=store+(row*(row+1))/2;
  246. }
  247.  
  248. void LowerTriangularMatrix::GetCol(MatrixRowCol& mrc)
  249. {
  250.    REPORT
  251.    int col=mrc.rowcol; mrc.skip=col; mrc.cw+=IsACopy; mrc.length=nrows; 
  252.    int i=nrows-col; mrc.storage=i; Real* ColCopy;
  253.    if ( !(mrc.cw*StoreHere) )
  254.    {
  255.       REPORT                                            // not accessed
  256.       ColCopy = new Real [i]; MatrixErrorNoSpace(ColCopy);
  257.       MONITOR_REAL_NEW("Make (LT GetCol)",i,ColCopy)
  258.       mrc.store = ColCopy-col;
  259.    }
  260.    else { REPORT ColCopy = mrc.store+col; }
  261.    if (+(mrc.cw*LoadOnEntry))
  262.    {
  263.       REPORT
  264.       Real* Mstore = store+(col*(col+3))/2;
  265.       while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
  266.    }
  267. }
  268.  
  269. void LowerTriangularMatrix::RestoreCol(MatrixRowCol& mrc)
  270. {
  271. //  if (mrc.cw*StoreOnExit)
  272.     {
  273.         REPORT
  274.         int col=mrc.rowcol; Real* Cstore = mrc.store+col;
  275.         Real* Mstore = store+(col*(col+3))/2; int i=nrows-col;
  276.       while (i--) { *Mstore = *Cstore++; Mstore += ++col; }
  277.    }
  278. }
  279.  
  280. void LowerTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrLT(); }
  281.                                      //712
  282. // routines for symmetric matrix
  283.  
  284. void SymmetricMatrix::GetRow(MatrixRowCol& mrc)
  285. {
  286.    REPORT                                                //571
  287.    mrc.skip=0; int row=mrc.rowcol; mrc.length=ncols;
  288.    if (+(mrc.cw*DirectPart))
  289.    {
  290.       REPORT
  291.       mrc.cw-=IsACopy; mrc.storage=row+1; mrc.store=store+(row*(row+1))/2;
  292.    }
  293.    else
  294.     {
  295.         mrc.cw+=IsACopy; mrc.storage=ncols;
  296.       Real* RowCopy = new Real [ncols]; MatrixErrorNoSpace(RowCopy);
  297.       MONITOR_REAL_NEW("Make (SymGetRow)",ncols,RowCopy) 
  298.       mrc.store = RowCopy;
  299.       if (+(mrc.cw*LoadOnEntry))
  300.       {
  301.             REPORT                                         // 544
  302.          Real* Mstore = store+(row*(row+1))/2; int i = row;
  303.          while (i--) *RowCopy++ = *Mstore++;
  304.          i = ncols-row;
  305.             while (i--) { *RowCopy++ = *Mstore; Mstore += ++row; }
  306.       }
  307.    }
  308. }
  309.  
  310. // need to check this out under StoreHere
  311.  
  312. void SymmetricMatrix::GetCol(MatrixRowCol& mrc)
  313. {
  314.     int col=mrc.rowcol; mrc.length=nrows;
  315.     if (mrc.cw >= StoreHere+DirectPart)
  316.     {
  317.         REPORT
  318.         mrc.skip=col; mrc.cw+=IsACopy;
  319.         int i=nrows-col; mrc.storage=i;
  320.         Real* ColCopy = mrc.store+col;
  321.         if (+(mrc.cw*LoadOnEntry))
  322.         {
  323.             REPORT                           // not accessed
  324.             Real* Mstore = store+(col*(col+3))/2;
  325.             while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
  326.         }
  327.     }
  328.     else
  329.     {
  330.         REPORT
  331.         mrc.skip=0;
  332.         if (+(mrc.cw*DirectPart))
  333.         {
  334.             REPORT
  335.             mrc.cw-=IsACopy; mrc.storage=col+1; mrc.store=store+(col*(col+1))/2;
  336.         }
  337.         else
  338.         {
  339.             mrc.cw+=IsACopy; mrc.storage=ncols; Real* ColCopy;
  340.             if ( !(mrc.cw*StoreHere) )
  341.             {
  342.                 REPORT                                      // not accessed
  343.                 ColCopy = new Real [ncols]; MatrixErrorNoSpace(ColCopy);
  344.                 MONITOR_REAL_NEW("Make (SymGetCol)",ncols,ColCopy)
  345.                 mrc.store = ColCopy;
  346.             }
  347.             else { REPORT ColCopy = mrc.store; }
  348.             if (+(mrc.cw*LoadOnEntry))
  349.             {
  350.                 REPORT
  351.                 Real* Mstore = store+(col*(col+1))/2; int i = col;
  352.                 while (i--) *ColCopy++ = *Mstore++;
  353.                 i = ncols-col;
  354.                 while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
  355.             }
  356.         }
  357.     }
  358. }
  359.  
  360. //void SymmetricMatrix::RestoreRow(int row, Real* Rstore)
  361. //{
  362. ////   if (cw*IsACopy && cw*StoreOnExit)
  363. //   {
  364. //      Real* Mstore = store+(row*(row+1))/2; int i = row+1;
  365. //      while (i--) *Mstore++ = *Rstore++;
  366. //   }
  367. //}
  368.  
  369. void SymmetricMatrix::RestoreCol(MatrixRowCol& mrc)
  370. {
  371. if (+(mrc.cw*StoreHere))
  372.     {
  373.         REPORT
  374.         // Really do restore column
  375.         // Restoring row instead would be possible but would make
  376.         // solve less efficient
  377.         // Same as LowerTriangular::RestoreCol
  378.  
  379.         int col=mrc.rowcol; Real* Cstore = mrc.store+col;
  380.         Real* Mstore = store+(col*(col+3))/2; int i = nrows-col;
  381.         while (i--) { *Mstore = *Cstore++; Mstore+= ++col; }
  382.     }
  383. }
  384.  
  385. // routines for row vector
  386.  
  387. void RowVector::GetCol(MatrixRowCol& mrc)
  388. {
  389.     REPORT
  390.     mrc.skip=0; mrc.storage=1; mrc.length=nrows;
  391.     if (mrc.cw >= StoreHere)
  392.    {
  393.       if (mrc.cw >= LoadOnEntry) { REPORT *(mrc.store) = *(store+mrc.rowcol); }
  394.       mrc.cw+=IsACopy;
  395.    }
  396.    else  { REPORT mrc.store = store+mrc.rowcol; mrc.cw-=IsACopy; }
  397.                                                          // not accessed
  398. }
  399.  
  400. void RowVector::NextCol(MatrixRowCol& mrc) 
  401. {
  402.    REPORT
  403.    if (+(mrc.cw*StoreHere))
  404.    {
  405.       if (+(mrc.cw*StoreOnExit)) { REPORT *(store+mrc.rowcol)=*(mrc.store); }
  406.                              // not accessed
  407.       mrc.rowcol++;
  408.       if (mrc.rowcol < ncols)
  409.       {
  410.      if (+(mrc.cw*LoadOnEntry)) { REPORT *(mrc.store)=*(store+mrc.rowcol); }
  411.       }
  412.       else { REPORT mrc.cw -= StoreOnExit; }
  413.    }
  414.    else  { REPORT mrc.rowcol++; mrc.store++; }             // not accessed
  415. }
  416.  
  417. void RowVector::RestoreCol(MatrixRowCol& mrc)
  418. {
  419.    REPORT                                            // not accessed
  420.    if (mrc.cw>=IsACopy)  { REPORT *(store+mrc.rowcol)=*(mrc.store); }
  421. }
  422.  
  423.  
  424. // routines for band matrices
  425.  
  426. void BandMatrix::GetRow(MatrixRowCol& mrc)
  427. {
  428.    REPORT
  429.    mrc.cw -= IsACopy; int r = mrc.rowcol; int w = lower+1+upper;
  430.    mrc.length=ncols;
  431.    int s = r-lower; mrc.store = store+(r*w-s); if (s<0) { w += s; s = 0; }
  432.    mrc.skip = s; s += w-ncols; if (s>0) w -= s; mrc.storage = w;
  433. }
  434.  
  435. // make special versions of this for upper and lower band matrices
  436.  
  437. void BandMatrix::NextRow(MatrixRowCol& mrc)
  438. {
  439.    REPORT
  440.    int r = ++mrc.rowcol; mrc.store += lower+upper;
  441.    if (r<=lower) mrc.storage++; else mrc.skip++;
  442.    if (r>=ncols-upper) mrc.storage--;
  443. }
  444.  
  445. void BandMatrix::GetCol(MatrixRowCol& mrc)
  446. {
  447.    REPORT
  448.    mrc.cw += IsACopy; int c = mrc.rowcol; int n = lower+upper; int w = n+1;
  449.    mrc.length=nrows;
  450.    int b; int s = c-upper; Real* ColCopy;
  451.    if (s<=0) { w += s; s = 0; b = c+lower; } else b = s*w+n;
  452.    mrc.skip = s; s += w-nrows; if (s>0) w -= s; mrc.storage = w;
  453.    if ( !(mrc.cw*StoreHere) )
  454.    {
  455.       REPORT
  456.       ColCopy = new Real [w]; MatrixErrorNoSpace(ColCopy);
  457.       MONITOR_REAL_NEW("Make (BMGetCol)",w,ColCopy)
  458.       mrc.store = ColCopy-mrc.skip;
  459.    }
  460.    else { REPORT ColCopy = mrc.store+mrc.skip; }
  461.    if (+(mrc.cw*LoadOnEntry))
  462.    {
  463.       REPORT
  464.       Real* Mstore = store+b;
  465.       while (w--) { *ColCopy++ = *Mstore; Mstore+=n; }
  466.    }
  467. }
  468.  
  469. void BandMatrix::RestoreCol(MatrixRowCol& mrc)
  470. {
  471. //  if (mrc.cw*StoreOnExit)
  472.     REPORT
  473.     int c = mrc.rowcol; int n = lower+upper; int s = c-upper;
  474.     Real* Mstore = store + ((s<=0) ? c+lower : s*n+s+n);
  475.     Real* Cstore = mrc.store+mrc.skip; int w = mrc.storage;
  476.     while (w--) { *Mstore = *Cstore++; Mstore += n; }
  477. }
  478.  
  479. // routines for symmetric band matrix
  480.  
  481. void SymmetricBandMatrix::GetRow(MatrixRowCol& mrc)
  482. {
  483.    REPORT
  484.    int r=mrc.rowcol; int s = r-lower; int w1 = lower+1; int o = r*w1;
  485.    mrc.length = ncols;
  486.    if (s<0) { w1 += s; o -= s; s = 0; }  mrc.skip = s;
  487.  
  488.    if (+(mrc.cw*DirectPart))
  489.       { REPORT  mrc.cw -= IsACopy; mrc.store = store+o-s; mrc.storage = w1; }
  490.    else
  491.    {
  492.       mrc.cw += IsACopy; int w = w1+lower; s += w-ncols;
  493.       if (s>0) w -= s; mrc.storage = w; int w2 = w-w1;
  494.       Real* RowCopy = new Real [w]; MatrixErrorNoSpace(RowCopy);
  495.       MONITOR_REAL_NEW("Make (SmBGetRow)",w,RowCopy) 
  496.       mrc.store = RowCopy-mrc.skip;
  497.       if (+(mrc.cw*LoadOnEntry))
  498.       {
  499.      REPORT
  500.          Real* Mstore = store+o;
  501.          while (w1--) *RowCopy++ = *Mstore++;   Mstore--;
  502.          while (w2--) { Mstore += lower; *RowCopy++ = *Mstore; }
  503.       }
  504.    }
  505. }
  506.  
  507. // need to check this out under StoreHere
  508.  
  509. void SymmetricBandMatrix::GetCol(MatrixRowCol& mrc)
  510. {
  511.    int c=mrc.rowcol; int w1 = lower+1; mrc.length=nrows;
  512.    if (mrc.cw >= StoreHere+DirectPart)
  513.    {
  514.       REPORT
  515.       mrc.cw += IsACopy;
  516.       Real* ColCopy; int b = c*w1+lower;
  517.       mrc.skip = c; c += w1-nrows; w1 -= c; mrc.storage = w1;
  518.       REPORT ColCopy = mrc.store+mrc.skip;
  519.       if (+(mrc.cw*LoadOnEntry))
  520.       {
  521.          REPORT
  522.          Real* Mstore = store+b;
  523.          while (w1--) { *ColCopy++ = *Mstore; Mstore += lower; }
  524.       }
  525.    }
  526.    else
  527.    {
  528.       REPORT
  529.       int s = c-lower; int o = c*w1;
  530.       if (s<0) { w1 += s; o -= s; s = 0; }  mrc.skip = s;
  531.  
  532.       if (+(mrc.cw*DirectPart))
  533.          { REPORT  mrc.cw -= IsACopy; mrc.store = store+o-s; mrc.storage = w1; }
  534.       else
  535.       {
  536.          mrc.cw += IsACopy; int w = w1+lower; s += w-ncols;
  537.          if (s>0) w -= s; mrc.storage = w; int w2 = w-w1; Real* ColCopy;
  538.          if ( !(mrc.cw*StoreHere) )
  539.          {
  540.             REPORT ColCopy = new Real [w]; MatrixErrorNoSpace(ColCopy);
  541.             MONITOR_REAL_NEW("Make (SmBGetCol)",w,ColCopy)
  542.             mrc.store = ColCopy-mrc.skip;
  543.          }
  544.          else { REPORT ColCopy = mrc.store+mrc.skip; }
  545.          if (+(mrc.cw*LoadOnEntry))
  546.          {
  547.             REPORT
  548.             Real* Mstore = store+o;
  549.             while (w1--) *ColCopy++ = *Mstore++;   Mstore--;
  550.             while (w2--) { Mstore += lower; *ColCopy++ = *Mstore; }
  551.          }
  552.       }
  553.    }
  554. }
  555.  
  556. void SymmetricBandMatrix::RestoreCol(MatrixRowCol& mrc)
  557. {
  558. if (+(mrc.cw*StoreHere))
  559.    {
  560.       REPORT
  561.       int c = mrc.rowcol;
  562.       Real* Mstore = store + c*lower+c+lower;
  563.       Real* Cstore = mrc.store+mrc.skip; int w = mrc.storage;
  564.       while (w--) { *Mstore = *Cstore++; Mstore += lower; }
  565.    }
  566. }
  567.  
  568.  
  569.